%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 蚁群算法
% 简介：蚁群算法寻找最大值
% 作者：Zhaojiang
% 日期：2023/10/10
% 企鹅：277746470
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% target fucntion
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
pause(2);
close all;clear;clc;
%% 设定算法参数
rng(1234);
var_num = 2;% 2个求解变量
lower_bound = -4;% 解的最小值（小数后4位）
upper_bound = 4;% 解的最大值（小数后4位）
bound_length = upper_bound-lower_bound;
ant_num = 100;% 蚂蚁数量
max_iteration = 200;% 最大迭代次数
Rho = 0.9;% 信息素挥发因子
P0 = 0.5;% 转移概率常数
step_size = 0.01;% 优化步长
global_amplitude = 0.5;
%% 算法数据初始化
ant = lower_bound +bound_length*rand(ant_num,var_num);% 随机生成蚁群
Tau = target(ant(:,1),ant(:,2));% 蚁群信息素初始化
[best_func_value,indexI] = max(Tau);% 最大蚁群信息素
best_var = ant(indexI,:);% 最佳变量值
best_Tau = best_func_value;% 最佳函数值
all_best_func_value =zeros(1,max_iteration);% 记录最佳函数值变化
%% 开始迭代求解
for iteration=1:max_iteration
    lamba = 1/iteration;% 启发概率(路径倒数)
    for i = 1:ant_num
        P = (best_Tau-Tau(i))/best_Tau;% 状态转移概率生成新解
        if P < P0% 局部搜索
            new_var =ant(i,:) + (2*rand(1,2)-1)*step_size*lamba;
        else% 全局搜索
            new_var = ant(i,:) + (2*rand(1,2)-1)*bound_length*global_amplitude;
        end
        % 变量限幅
        for j = 1:var_num
            if new_var(j)>upper_bound
                new_var(j) = upper_bound;
            elseif new_var(j)<lower_bound
                 new_var(j) = lower_bound;           
            end
        end
        % 新解函数值
        func_value = target(new_var(1),new_var(2));
        if func_value > Tau(i)% 新解有效
            ant(i,:) = new_var;
        end
        % 是否保存最佳值信息
        func_value = target(ant(i,1),ant(i,2));
        if func_value >best_func_value
            best_var = ant(i,:);
            best_func_value = func_value;
            best_iteration = iteration;
        end
        % 信息素挥发
        Tau(i) = (1 - Rho) * Tau(i) + target(ant(i,1),ant(i,2));
        % 最佳信息素更新
        if best_Tau < Tau(i)
            best_Tau = Tau(i);
        end
    end 
    % 最佳值记录
    all_best_func_value(iteration) = best_func_value;
end
%% 处理结果
% 打印结果
max_x = best_var(1);
max_y = best_var(2);
max_ans = best_func_value;
disp 最优迭代次数
disp(best_iteration)
disp 最大值位置
fprintf("x=%f\ny=%f\nz=%f\n",max_x,max_y,max_ans);
% 可视化结果
figure('position',[200 200 1280 480])
subplot(1,2,1);hold on;
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
view(-37.5,30)
pause(0.1)
scatter3(max_x,max_y,max_ans,'r*');

subplot(1,2,2);hold on;
x = 1:max_iteration;
y1= all_best_func_value;
try
    assert(length(x)==length(y1))
    plot(x,y1,'-','LineWidth',1.2)
    xlabel('x');ylabel('y');
    title('迭代数据')
catch
    disp("plot failed")
    disp(size(x))
    disp(size(y1))
end